Combining measurements from different sensors

ABSTRACT

A method for combining measurements from two or more independent measurement channels, particularly physiological measurements such as heart rate. Independent measurements of heart rate, for instance by ECG and pulse oximetry, can be combined to derive an improved measurement eliminating artefacts on one channel. A model of the process generating the physiological parameter, e.g., the heart rate, is constructed and is run independently for each channel to generate predictions of the parameter. The measured values are compared with the predicted values and the differences are used as an indication of the confidence in the measurement. The measurements from the two channels are ombined using weights calculated from the respective differences.

This application is the U.S. national phase of international applicationPCT/GB01/02544 filed Jun. 8, 2001, which designated the U.S.

BACKGROUND AND SUMMARY

The present invention relates to a method and apparatus for combiningmeasurements from different sensors in order to provide an improvedmeasurement of a parameter. It is particularly applicable to themeasurement of physiological parameters.

Certain parameters can be measured in more than one way. This is usefulin giving independent measures of the same quantity. For instance, inthe medical field the heart rate can be measured both from anelectrocardiogram (ECG) and from a pulse oximetry waveform (used tocalculate oxygen saturation). drawings FIG. 4 illustrates schematicallythese-two-waveforms, FIG. 4 a being the electrocardiogram with the heartrate illustrated as HR₁, and FIG. 4 b the pulse oximetry waveform withthe heart rate illustrated as HR₂. The heart rate is a parameter whichcan undergo sudden changes. Some of these changes are validphysiological changes, for example ectopic beats which occurprematurely, and therefore give rise to a temporary increase in theheart rate. FIG. 5 illustrates the occurrence of an ectopic beat 50found in both the electrocardiogram trace and the pulse oximetrywaveform. The shorter interval between the preceeding beat and theectopic beat 50 manifests itself in a measurement of the heart rate as asudden increase in the heart rate. FIGS. 1 and 2 of the accompanyingdrawings show time plots of the heart rate measured by pulse oximetry(FIG. 1) and ECG (FIG. 2). It can be seen that in FIGS. 1 and 2 theheart rate in. the early part of the plot is generally of the order of80 bpm, but that there are occasional sudden increases in heart rate,such as indicated at 10 and 20 which are caused by ectopic beats andthus appear both in the measurement by pulse oximetry and themeasurement by ECG.

However, in addition to changes in the measured heart rate deriving fromvalid physiological changes, other changes occur which are notphysiologically valid, for instance being caused by sudden movement ofthe sensors on the body surface (e.g. chest movement with ECGelectrodes). FIG. 6 illustrates the presence of artefacts 60 on thepulse oximetry waveform which shorten the interval between apparentbeats and thus result in apparent increases in the heart rate. Thesechanges are reflected in one measurement, but not the other, asindicated at 12 and 22 in FIGS. 1 and 2 respectively. The fact that thechanges appear in one measurement but not the other means that the twomeasurements could be combined to help decide which heart rate changesare valid physiological ones, and which are artefacts. However, thenormal approach of validating one measurement channel against the otherinvolving cross-correlation of the two measurements invariably failsbecause it is not possible to know in advance (for each recording, foreach patient) what value to give to the threshold for accepting, ratherthan rejecting a change in the heart rate as being valid. Thus althoughit would appear from FIGS. 1 and 2 that a threshold could be set whichwould eliminate changes such as indicated as 22, such a threshold is notappropriate for all patients for all recordings, and does not help withthe pulse oximetry waveform. The problems are increased in the event ofatrial fibrillation when the heart rate changes rapidly as indicated inthe region AF in FIGS. 1, 2 and 3.

Similar problems arise in other fields where a parameter is measured viatwo or more measurement channels.

The present invention provides a method and apparatus for improvingmeasurement of a parameter by combining two measurements of it in a waywhich allows valid changes to be distinguished from artefacts.Accordingly it provides a method of measuring a parameter comprising thesteps of: predicting the value of each of two measurements of theparameter, making two measurements of the parameter to produce twomeasured values of the parameter, calculating the respective differencesbetween the predicted values and the measured values, and combining thetwo measured values with weights determined by said differences.

Thus with the present invention a prediction is made for eachmeasurement and the actual measurement is compared with its prediction.The difference is computed, which is termed the “innovation”, and thisinnovation is used to calculate a weight which will be given to thatmeasurement when it is combined with the other measurement, alsoweighted according to its innovation. The weights are calculated so thatif the innovation on one measurement channel is high, whereas theinnovation on the other measurement channel is low, the measurement fromthe low innovation channel is more heavily weighted. This is because ahigh level of innovation from one channel coinciding with a lowinnovation on the other channel is regarded as indicative of an artefacton the higher innovation channel. Thus, the weight given to each valuewhen the values are combined is inversely related to the square ormodulus of the difference between the measured value and its predictedvalue.

In one embodiment the measured values can be combined according to theformula:

$\begin{matrix}{M = {{M_{1}\frac{\sigma_{2}^{2}}{\sigma_{1}^{2} + \sigma_{2}^{2}}} + {M_{2}\frac{\sigma_{1}^{2}}{\sigma_{1}^{2} + \sigma_{2}^{2}}}}} & (1)\end{matrix}$where M₁ and M₂ are the two measured values, and σ₁ and σ₂ are thedifferences between the two measured values and their respectivepredicted values.

The steps of prediction, measurement, calculation and combination arepreferably repeated continuously, with the predicted value for each ofthe measurements being based on a linear predictive model, e.g thepredicted value is based on its preceding predicted value and thepreceding innovation (i.e. the difference between the precedingpredicted value and the preceding measurement). The predicted value canbe obtained by adding to the preceding predicted value a constant timesthe innovation. The constant is preferably a positive value less than orequal to unity. Alternatively the predicted value for each of the twoindependent measurements can be calculated by using. a non-linear,predictive model such as a neural network.

In one embodiment the predicted values can be based on a mathematicalmodel of the system, which may include estimates for process noise andsensor (measurement) noise. Two independent models may be used, one foreach of the measurement channels, and the models can include estimatesfor the process noise and sensor noise, which can be the same for thetwo channels. In one embodiment the models are Kalman filters.

The method is particularly applicable to the measurement of heart rate,in which case the two measurement channels can be from anelectrocardiogram and a pulse oximetry waveform, though it is applicableto any other measurement of a parameter which can be derived from two ormore sources. Thus the method is applicable for more than twomeasurement channels, and both where the measurements are independentand where they are not truly independent such as from multiple leads ofan ECG.

The invention can also provide for detection of movement artefacts. Inthis instance high values of innovation are obtained on both channelsfor the period of movement, and this can be used as a trigger to discardthe sections of data which are corrupted by that movement.

It will be appreciated that the invention can be embodied using computersoftware and thus the invention extends to a computer program forcontrolling and executing the method or parts of it, and to a computerreadable storage medium carrying the program. The invention also extendsto corresponding apparatus for carrying out the method.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention will be further described by way of non-limitative examplewith reference to the accompanying drawings in which:

FIG. 1 illustrates a plot of heart rate measured by pulse oximetry;

FIG. 2 illustrates a plot of heart rate measured by an ECG;

FIG. 3 illustrates the result of combining the two plots of FIGS. 1 and2 according to an embodiment of the invention;

FIG. 4 illustrates schematically heart beats on an ECG and pulseoximetry waveform;

FIG. 5 illustrates schematically ectopic beats appearing on an ECG andpulse oximetry waveform;

FIG. 6 illustrates an ECG trace and pulse oximetry waveform withartefacts on the pulse oximetry waveform;

FIG. 7 illustrates a plot of heart rate measured by an ECG;

FIG. 8 illustrates predicted values for the heart rate according to oneembodiment of the invention;

FIG. 9 illustrates the innovation obtained from FIGS. 7 and 8;

FIG. 10 illustrates the variance obtained from FIGS. 7 and 8.

DETAILED DESCRIPTION OF EXAMPLE EMBODIMENTS

An embodiment of the invention will now be described in which theinvention is applied in the medical field to the measurement of heartrate using ECG and pulse oximetry. As illustrated in FIG. 4a the heartrate measured by ECG is derived from the interval between two successiveR-wave peaks. The heart rate measurement derived from the pulse oximetrywaveform is obtained from the interval between two successive peaks (ortroughs) as illustrated in FIG. 4 b. FIGS. 1 and 2 illustrate heart rateplots derived from these two measurements.

With this embodiment of the present invention a model of the processgenerating the heart rate is constructed. The same model is runindependently for each measurement source (i.e. one for the ECGmeasurement channel and one for the pulse oximetry measurement channel).In this embodiment the model is a Kalman filter. In general a Kalmanfilter uses a process model and an observation model. The process modelmodels the state of the system at time t+1 in terms of its state at timet. The measurement or observation model indicates how the measurement attime t is related to the state of the system at time t. Thus in generalterms:

-   -   X_(t+1)=Ax₁+w (Process model)    -   Y_(t)=Cx_(t)+v (Observation model)        where:    -   w˜N (O, Q)— Gaussian process noise with zero mean and variance Q    -   v˜N (O, R)— Gaussian measurement noise with zero mean and        variance R        and:    -   k—vector of state variables x    -   n—vector of observable or measurements y    -   State x evolves according to simple first-order Markov dynamics;    -   A is the k×k state transition matrix    -   Each measurement vector y is related to the current state by a        linear observation process; C is the n×k observation or        measurement matrix

In this embodiment the general Kalman filter is simplified to use scalarquantities and the same process and measurement noise models (w, v) areused on both measurement channels. Thus the simplified Kalman filter isas follows:

-   -   x_(t+1)=Ax_(t)+w (Process model)    -   y_(t)=Cx_(t)+v (Observation model)

The model is further simplified by setting C=1, implying that the heartrate is both the state describing the process and the measurement.Further, it is assumed that A=1, implying that the next heart rate isthe same as the previous one with the variability allowed for by theprocess noise model.

Using this model, on each channel, the process of combining the twomeasurements then involves the following steps:

-   -   1. From knowledge of previous history up to time t (a) predict        the next state xpred; (b) from xpred predict the next        measurement ypred    -   2. Make the measurement y_(t+1)    -   3. Compute the innovation:        -   y_(t+1)=ypred+ε_(t+1)        -   where ε_(t+1) is the difference between the actual value and            the predicted value: the innovation.    -   4. Compute the variance σ_(t+1) ²:        σ_(t+1) ²=ε_(t+1) ²        -   σ² the variance, is the inverse of the “confidence” which is            associated with the prediction

-   5. Mix the heart rate measurements in inverse proportion to the    variance associated with each one:

${HR} = {{{HR}_{1}\frac{\sigma_{2}^{2}}{\sigma_{1}^{2} + \sigma_{2}^{2}}} + {{HR}_{2}\frac{\sigma_{1}^{2}}{\sigma_{1}^{2} + \sigma_{2}^{2}}}}$An example of an implementation of this model in MATLAB is given inAppendix 1. That example is general, and will work for vectorquantities, though in this embodiment the quantities are scalar. It canbe seen from appendix 1 that the predicted heart rate for each newmeasurement cycle (xnew) is equal to the previously predicted value(xpred) plus the Kalman gain K times the innovation e. The Kalman gain Kis derived from the predicted variance Vpred and the measurementvariance R. The predicted variance is derived from the previouspredicted variance and the process noise Q. To start the process off itis initialized using an initial value of the heart rate as 80 and aninitial value of the state variance of 100. The process noise in the Qin this embodiment is set to 5 and the measurement noise variance R isset to 10.

It will be clear from the implementation that, as normal with a Kalmanfilter, the variance and Kalman gain are not dependent on themeasurement values. The measurement values are only used in the newprediction of heart rate via the innovation e. Thus it will be notedthat for the constant values of Q and R used in this example the Kalmangain K tends to 0.5 and the state co-variance V tends to 5. However, Kcan be made adaptive by modifying the values for the variance constantsQ and R, preferably the process variance Q, according to the type ofprocess being encountered, for example atrial fibrillation (where thereis a high level of process noise) as opposed to a healthy heart rate(during which there is a low level of process noise).

FIGS. 8, 9 and 10 illustrate values for the estimated heart rate xpred,the innovation e and the variance σ² for the heart rate plot shown inFIG. 7 (in this case the ECG measurement channel).

FIG. 3 illustrates the result of combining the two measurement plotsfrom FIGS. 1 and 2 using this embodiment. It can be seen that themovement artefacts 22 on the ECG channel in FIG. 2 have been removedfrom the combined measurement, even though they occur during a period ofatrial fibrillation.

Thus with this embodiment the difference between the measurement and itspredicted value is used to indicated the degree of confidence in thatmeasurement. The higher the difference the lower the confidence. Formula(1) above-is used to combine the two measurements. This can besummarised as follows:

-   -   1) Valid heart rates on both channels: low innovation on both        channels; weight both measurements equally    -   2) Valid sudden change (e.g. ectopic beat) seen on both        channels: high innovation but on both channels; therefore,        measurements are again weighted equally    -   3) Artefact on one channel: high innovation on one channel only;        therefore the information from that channel is ignored by being        given a very low weighting (low confidence).

The method can also be used to provide a movement artefact detector,i.e. to detect when movement artefact is present on both channels andtherefore no useful information is available. This is characterised byhigh values of innovation on both channels for a sustained period oftime. This can be used to discard the sections of raw data which arecorrupted by this movement and to indicate that no valid heart rateestimate can be derived during those periods.

Appendix 1 load -ascii ecghr_13 data_file = ecghr_13; time =data_file(;,1); hr = data_file(:,2); start = 1 stop =size(data_file),1); fprintf(‘number of R—R intervals detected = %d \n’,stop); hr_limit = 200; % X(t + 1) = A X(t) + noise(Q) - process modelwith Q as variance of noise w % Y(t) = C X(t) + noise(R) - measurementmodel with R as variance of noise v ss = 1; % state size - sets to onedimensional, ie scalar though routine works for vectors os = 1; %observation size - sets to one dimensional, ie scalar A = [1]; % assumex(t + 1) = x(t) C = [1]; % assume y = x Q = 5.0*eye(ss); % processnoise - eye is the identity matrix in MATLAB-here just unity R =10.0*eye(os); % measurement noise variance initx = [80]; % initial statevalue (HR of 80 bpm) initV = 100*eye(ss); % initial state variance xnew= initx; % - initialisation Vnew = initV; for i = start:stop %-start ofcycle x = xnew; % update from previous cycle V = Vnew; % update fromprevious cycle xpred = A*x; % prediction of state Vpred = A*V*A′ + Q; %prediction of state covariance, A′ is transpose of A ypred(i) = C*xpred;% prediction of measurement y(i) = hr(i); % “make measurement” e = y(i)− ypred(i); % calculate innovation innov(i) = e; % for plottingsigma2(i) = e*e; % variance for saving S = C*Vpred*C′ + R; % innovationcovariance Sinv = inv(S); % invert S K = Vpred*C′ *Sinv; % computeKalman gain xnew = xpred + K*e; % update state by the innovationcontrolled by the Kalman gain Vnew = (eye (ss) − K*C) *Vpred; % updatestate covariance end end

1. A method of measuring a parameter comprising using a suitablyprogrammed computer to perform steps including: predicting a value ofeach of two independent measurements of the parameter for twoindependent measurement channels, making two independent measurements ofthe parameter via said two independent measurement channels to producetwo measured values of the parameter, calculating respective differencesbetween the predicted values and the measured values for each of saidtwo independent measurement channels, combining the two measured valueswith weights which vary inversely with said respective differences toproduce a combined value, and outputting the combined value to a user.2. A method according to claim 1 in which the predicting, the making oftwo independent measurements, the calculating and the combining arerepeated, the predicted value for each of the two independentmeasurements being based on a preceding predicted value and a differencebetween the preceding predicted value and a preceding measured value. 3.A method according to claim 2 in which the predicted value for each ofthe two independent measurements is calculated by using a linearpredictive model for each of said two independent measurement channels.4. A method according to claim 2 in which the predicted value for eachof the two independent measurements is calculated by using a non-linearpredictive model for each of said two independent measurement channels.5. A method according to claim 3 in which each model is adaptive, and itadapts in dependence upon the amount of process noise in the respectivemeasurement channels.
 6. A method according to claim 1 in which in thestep of combining the two measured values the weights are inverselyproportional to the squares of said respective differences.
 7. A methodof measuring a parameter comprising using a suitably programmed computerto perform steps including: predicting a value of each of twoindependent measurements of the parameter for two independentmeasurement channels, making two independent measurements of theparameter via said two independent measurement channels to produce twomeasured values of the parameter, calculating respective differencesbetween the predicted values and the measured values for each of saidtwo independent measurement channels, combining the two measured valueswith weights determined by said respective differences to produce acombined value M, and outputting the combined value to a user, whereinthe two measured values are combined according to the formula:—$M = {{M_{1}\frac{\sigma_{2}^{2}}{\sigma_{1}^{2} + \sigma_{2}^{2}}} + {M_{2}\frac{\sigma_{1}^{2}}{\sigma_{1}^{2} + \sigma_{2}^{2}}}}$where M1 and M2 are the two measured values, and σ₁ and σ₂ aredifferences between the two measured values and their respectivepredicted values.
 8. A method according to claim 1 in which thepredicted values for the respective measurements are based on respectivemodels for measurement of said parameter via said two independentmeasurement channels.
 9. A method according to claim 8 in which themodels include estimates for process noise and sensor noise for said twoindependent measurement channels.
 10. A method according to claim 8 inwhich the respective models are mutually independent.
 11. A methodaccording to claim 10 in which the respective models include the sameestimates for process noise and sensor noise.
 12. A method of measuringa parameter comprising using a suitably programmed computer to performsteps including: predicting a value of each of two independentmeasurements of the parameter for two independent measurement channelsusing respective models for measurement of the parameter, making twoindependent measurements of the parameter via said two independentmeasurement channels to produce two measured values of the parameter,calculating respective differences between the predicted values and themeasured values for each of said two independent measurement channels,combining the two measured values with weights which vary inversely withsaid respective differences to produce a combined value, and outputtingthe combined value to a user, wherein the respective models compriseKalman filters.
 13. A method according to claim 1 further comprisingdiscarding series of measurements for which the differences between bothmeasured values and their predicted values exceed a predeterminedthreshold for a predetermined period of time.
 14. A method according toclaim 1 in which the parameter is heart rate.
 15. A method according toclaim 14 in which the two independent measurements are made from anelectrocardiograph and a pulse oximetry waveform respectivelyconstituting said two independent measurement channels.
 16. A methodaccording to claim 14 in which the two measurements are made from amultiple lead ECG recording.
 17. A method according to claim 1 in whichthere are more than two measurements.
 18. A method according to claim 1further comprising identifying movement artifacts based on the values ofthe differences between both measured values and their predicted values.19. A computer program encoded on a computer-readable storage medium andcomprising program code which, when executed, performs the method ofclaim
 1. 20. Apparatus configured to perform the method of claim
 1. 21.A method according to claim 4 in which each model is adaptive, and itadapts in dependence upon the amount of process noise in the respectivemeasurement channels.
 22. A method of measuring a physiologicalparameter comprising using a suitably programmed computer to performsteps including: predicting by use of respective Kalman filters a valueof each of two independent measurements of the physiological parameterfor two independent measurement channels, making two independentmeasurements of the physiological parameter via the two independentmeasurement channels to produce two measured values of the physiologicalparameter, calculating respective differences between the predictedvalues and the measured values for each of the two independentmeasurement channels, combining the two measured values with weightswhich vary inversely with which vary inversely with the respectivedifferences to produce a physiological parameter measurement, andoutputting the physiological parameter measurement to a user.
 23. Themethod according to claim 22, in which the physiological parameter isheart rate.
 24. The method according to claim 22, in which a first ofthe two measurement channels corresponds to an electrocardiograph and asecond of the two measurement channels corresponds to pulse oximetry.25. A computer-readable storage medium on which computer-executableinstructions for performing the method according to claim 22 areencoded.
 26. A method of measuring a parameter comprising using asuitably programmed computer to perform steps including: predicting avalue of each of two independent measurements of the parameter for twoindependent measurement channels; making two independent measurements ofthe parameter via said two independent measurement channels to producetwo measured current values of the parameter; calculating respectivedifferences between the predicted values and the measured current valuesfor each of said two independent measurement channels; and combining thetwo measured current values with weights which vary inversely with saidrespective differences to produce and output to a user a measurement ofthe parameter which is based on the two measured current values.
 27. Amethod of measuring a parameter comprising using a suitably programmedcomputer to perform steps including: predicting a value of each of twoindependent measurements of the parameter for two independentmeasurement channels; making two independent measurements of theparameter via said two independent measurement channels to produce twomeasured current values of the parameter; calculating respectivedifferences between the predicted values and the measured current valuesfor each of said two independent measurement channels; combining the twomeasured current values with weights which vary inversely with saidrespective differences to produce a measurement of the parameter whichis based on the two measured current values; and displaying themeasurement of the parameter.
 28. A method according to claim 27 whereinthe measurement of the parameter is displayed as part of a plot.
 29. Acomputer-readable storage medium on which computer-executableinstructions for performing the method according to claim 27 areencoded.